home *** CD-ROM | disk | FTP | other *** search
/ United Public Domain Gold 2 / United Public Domain Gold 2.iso / utilities / pu309.dms / pu309.adf / Sat_Tracking / sat-parta.c < prev    next >
C/C++ Source or Header  |  1990-09-29  |  14KB  |  543 lines
  1.  
  2. /* Change Log
  3.     5/10/86     v2.2 Added single character aliases for satellite
  4.             names.
  5.  
  6.     4/30/86     v2.1 Print blank line if satellite dips below
  7.             horizon and reappears during same orbit and day
  8.  
  9.     4/29/86     v2.0 Changed GetSatelliteParams() to use AMsAT's
  10.             "kepler.dat" file. Moved schedule to "mode.dat" file.
  11.  
  12.     4/22/86     v1.3  Inserted N8FJB's suggestions for variable naming
  13.             which maintain 8 character uniqueness.
  14.             Also removed "include" file orbitr.h, which had two
  15.             definitions of external functions defined in orbit.c
  16.                 -K3MC 
  17.  
  18.     4/1/86        v1.2  Corrected a scanf conversion to %d for an int
  19.             type.     -K3MC
  20.  
  21.     3/19/86     v1.1  Changed GetSatelliteParams to not pass NULL
  22.             to sscanf.
  23.  
  24.                                     */
  25.  
  26. #include "df0:include/stdio.h"
  27. #include "df0:include/math.h"
  28.  
  29. extern double Kepler();
  30. extern double GetDay();
  31.  
  32. #define PI 3.14159265
  33. #define PI2 (PI*2)
  34. #define MinutesPerDay (24*60.0)
  35. #define SecondsPerDay (60*MinutesPerDay)
  36. #define HalfSecond (0.5/SecondsPerDay)
  37. #define EarthRadius 6378.16        /* Kilometers        */
  38. #define C 2.997925e5            /* Kilometers/Second    */
  39. #define DegreesPerRadian (180/PI)
  40. #define RadiansPerDegree (PI/180)
  41. #define ABS(x) ((x) < 0 ? (-(x)) : (x))
  42. #define SQR(x) ((x)*(x))
  43.  
  44. #define MaxModes 10
  45. typedef struct {
  46.         int MinPhase,MaxPhase;
  47.         char ModeStr[20];
  48.            }  ModeRec;
  49.   
  50. static char Str1[30] ={ "N3EMO Orbit Simulator  v2.2"};
  51. static char Str2[30] = {"N3DPU Version for Amiga "};
  52.    
  53.  
  54.    /*    Keplerian Elements and misc. data for the satellite         */
  55.  
  56.     double  EpochDay;            /* time of epoch         */
  57.     double  EpochMeanAnomaly;    /* Mean Anomaly at epoch     */
  58.     long    EpochOrbitNum;         /* Integer orbit # of epoch     */
  59.     double  EpochRAAN;            /* RAAN at epoch         */
  60.     double  epochMeanMotion;        /* Revolutions/day         */
  61.     double  OrbitalDecay;        /* Revolutions/day^2         */
  62.     double  EpochArgPerigee;        /* argument of perigee at epoch  */
  63.     double  Eccentricity;
  64.     double  Inclination;
  65.     char    SatName[100];
  66.     int     ElementSet;
  67.     double  BeaconFreq;            /* Mhz, used for doppler calc     */
  68.     double  MaxPhase;            /* Phase units in 1 orbit     */
  69.     double  perigeePhase;
  70.     int     NumModes;
  71.     ModeRec Modes[MaxModes];
  72.     int     PrintApogee;
  73.  
  74.     /* Simulation Parameters */
  75.  
  76.     double StartTime,EndTime, StepTime; /* In Days, 1 = New Year    */
  77.                     /*    of reference year    */
  78.  
  79.     /* Site Parameters */
  80.    
  81.      char SiteName[100];
  82.     double SiteLat,SiteLong,SiteAltitude,SiteMinElev;
  83.  
  84.  
  85. /* List the satellites in kepler.dat, and return the number found */
  86. ListSatellites()
  87. {
  88.     char str[100];
  89.     FILE *InFile;
  90.     char satchar;
  91.     int NumSatellites;
  92.  
  93.     printf("Available satellites:\n");
  94.  
  95.     if ((InFile = fopen("kepler.dat","r")) == 0)
  96.     {
  97.     printf(" kepler.dat not found\n");    /* Bitwise Or  "|" */
  98.  
  99.     exit(-1);
  100.     }
  101.  
  102.     satchar = 'a';
  103.     NumSatellites = 0;
  104.     while (fgets(str,100,InFile))
  105.  if (sscanf(str,"Satellite: %s",str) == 1)
  106.         {
  107.         printf("    %c) %s\n",satchar++,str);
  108.         NumSatellites++;
  109.         }
  110.  
  111.     fclose(InFile);
  112.  
  113.     return NumSatellites;
  114. }
  115.  
  116. /* Match and skip over a string in the input file. Exits on failure. */
  117.  
  118. MatchStr(InFile,FileName,Target)
  119. FILE *InFile;
  120. char *FileName,*Target;
  121. {
  122.     char str[100];
  123.  
  124.     fgets(str,strlen(Target)+1,InFile);
  125.     if (strcmp(Target,str))
  126.        {
  127.        printf("%s: found \"%s\" while expecting \"%s\n\"",FileName,str,Target);
  128.        exit(-1);
  129.        }
  130. }
  131.  
  132. GetSatelliteParams()
  133. {
  134.     FILE *InFile;
  135.     char str[100];
  136.     int EpochYear;
  137.    /*    double EpochHour,EpochMinute,EpochSecond;  */
  138.     int found;
  139.     int NumSatellites;
  140.     char satchar;
  141.  
  142.     NumSatellites = ListSatellites();
  143.  
  144.     found = 0;
  145.  
  146.     while (!found)
  147.     {
  148.     printf("Letter or satellite name :");
  149.     gets(SatName);
  150.  
  151.     if ((InFile = fopen("kepler.dat","r")) == 0)
  152.         {
  153.         printf("kepler.dat not found\n");
  154.         exit(-1);
  155.         }
  156.  
  157.     if (strlen(SatName) == 1)
  158.         {            /* use single character label */
  159.         satchar = SatName[0];
  160.         if (satchar < 'a' || satchar > 'a'+NumSatellites-1)
  161.         {
  162.         printf("'%c' is out of range\n",satchar);
  163.         fclose(InFile);
  164.         continue;
  165.         }
  166.  
  167.         do    
  168.         do  /* find line beginning with "Satellite: " */
  169.             fgets(str,100,InFile);
  170.         while (sscanf(str,"Satellite: %s",SatName) != 1);
  171.          while (satchar-- > 'a');
  172.         found = 1;
  173.         }
  174.         
  175.      else 
  176.          {
  177.          while (!found)  /* use satellite name */
  178.         {
  179.         if (! fgets(str,100,InFile))
  180.             break;    /* EOF */
  181.  
  182.         if (sscanf(str,"Satellite: %s",str) == 1
  183.             && strcmp(SatName,str) == 0)
  184.             found = 1;
  185.         }
  186.  
  187.         if (!found)
  188.         {
  189.         printf("Satellite %s not found\n",SatName);
  190.         fclose(InFile);
  191.         }
  192.         }
  193.     }
  194.  
  195.     BeaconFreq = 146.0;  /* Default value */
  196.  
  197.     fgets(str,100,InFile);    /* Skip line */
  198.  
  199.     MatchStr(InFile,"kepler.dat","Epoch time:");
  200.     fgets(str,100,InFile);
  201.     sscanf(str,"%lf",&EpochDay);
  202.  
  203.     EpochYear = EpochDay / 1000.0;
  204.     EpochDay -= EpochYear*1000.0;
  205.     EpochDay = GetDay(EpochYear,1,EpochDay);
  206.  
  207.    
  208.     MatchStr(InFile,"kepler.dat","Element set:");
  209.     fgets(str,100,InFile);
  210.     sscanf(str,"%d",&ElementSet);
  211.  
  212.     MatchStr(InFile,"kepler.dat","Inclination:");
  213.     fgets(str,100,InFile);
  214.     sscanf(str,"%lf",&Inclination);
  215.     Inclination *= RadiansPerDegree;
  216.  
  217.     MatchStr(InFile,"kepler.dat","RA of node:");
  218.     fgets(str,100,InFile);
  219.     sscanf(str,"%lf",&EpochRAAN);
  220.     EpochRAAN *= RadiansPerDegree;
  221.  
  222.     MatchStr(InFile,"kepler.dat","Eccentricity:");
  223.     fgets(str,100,InFile);
  224.     sscanf(str,"%lf",&Eccentricity);
  225.  
  226.     MatchStr(InFile,"kepler.dat","Arg of perigee:");
  227.     fgets(str,100,InFile);
  228.     sscanf(str,"%lf",&EpochArgPerigee);
  229.     EpochArgPerigee *= RadiansPerDegree;
  230.  
  231.     MatchStr(InFile,"kepler.dat","Mean anomaly:");
  232.     fgets(str,100,InFile);
  233.     sscanf(str,"%lf",&EpochMeanAnomaly);
  234.     EpochMeanAnomaly *= RadiansPerDegree;
  235.  
  236.     MatchStr(InFile,"kepler.dat","Mean motion:");
  237.     fgets(str,100,InFile);
  238.     sscanf(str,"%lf",&epochMeanMotion);
  239.  
  240.     MatchStr(InFile,"kepler.dat","Decay rate:");
  241.     fgets(str,100,InFile);
  242.     sscanf(str,"%lf",&OrbitalDecay);
  243.  
  244.     MatchStr(InFile,"kepler.dat","Epoch rev:");
  245.     fgets(str,100,InFile);
  246.     sscanf(str,"%ld",&EpochOrbitNum);
  247.  
  248.     while (1)
  249.         {
  250.         if (! fgets(str,100,InFile))
  251.         break;    /* EOF */
  252.         if (strlen(str) <= 2)
  253.         break;    /* Blank line */
  254.         if (sscanf(str,"Beacon: %lf",&BeaconFreq) == 1)
  255.         break;    /* found it */
  256.         }
  257.  
  258.  
  259.     PrintApogee = (Eccentricity >= 0.3);
  260.  
  261.     perigeePhase = 0; MaxPhase = 360; /* Default values */
  262.     NumModes = 0;
  263.  
  264.     if ((InFile = fopen("mode.dat","r")) == 0)
  265.     return;
  266.  
  267.     found = 0;
  268.     while (!found)
  269.     {
  270.     if (! fgets(str,100,InFile))
  271.         break;    /* EOF */
  272.     if (sscanf(str,"Satellite: %s",str) == 1
  273.         && strcmp(SatName,str) == 0)
  274.         found = 1;
  275.     }
  276.     
  277.     if (found)
  278.     {
  279.     MatchStr(InFile,"mode.dat","Perigee phase:");
  280.     fgets(str,100,InFile);
  281.     sscanf(str,"%lf",&perigeePhase);
  282.  
  283.     MatchStr(InFile,"mode.dat","Max phase:");
  284.     fgets(str,100,InFile);
  285.     sscanf(str,"%lf",&MaxPhase);
  286.  
  287.  
  288.     while (fgets(str,100,InFile))
  289.         {
  290.         if (sscanf(str,"Mode: %20s from %d to %d",Modes[NumModes].ModeStr,
  291.         &Modes[NumModes].MinPhase,&Modes[NumModes].MaxPhase) == 3
  292.           && NumModes < MaxModes)
  293.           NumModes++;
  294.         }
  295.     fclose(InFile);
  296.     }
  297.  
  298.  
  299.  
  300. }
  301.  
  302.  
  303. GetSiteParams()
  304. {
  305.     FILE *InFile;
  306.     char str[100];
  307.  
  308. /*    printf("Site name :");   */
  309. /*    gets(str);    */
  310. /*    strcat(str,".sit");   */
  311.  
  312.     if ((InFile = fopen("Site.dat","r")) == 0)
  313.     {
  314.     printf("%s not found\n",str);
  315.     exit(-1);
  316.     }
  317.  
  318.     fgets(SiteName,100,InFile);
  319.  
  320.     fgets(str,100,InFile);
  321.     sscanf(str,"%lf",&SiteLat);
  322.     SiteLat *= RadiansPerDegree;
  323.  
  324.     fgets(str,100,InFile);
  325.     sscanf(str,"%lf",&SiteLong);
  326.     SiteLong *= RadiansPerDegree;
  327.  
  328.     fgets(str,100,InFile);
  329.     sscanf(str,"%lf",&SiteAltitude);
  330.     SiteAltitude /= 1000;   /* convert to km */
  331.  
  332.     fgets(str,100,InFile);
  333.     sscanf(str,"%lf",&SiteMinElev);
  334.     SiteMinElev *= RadiansPerDegree;
  335. }
  336.  
  337. GetSimulationParams()
  338. {
  339.     double hour,duration;
  340.     int Month,Day,Year;
  341.  
  342.     printf("Start date (UTC) (Month Day Year) :");
  343.     scanf("%d%d%d",&Month,&Day,&Year);
  344.  
  345.     StartTime = GetDay(Year,Month,(double) Day);
  346.  
  347.     printf("Starting Hour (UTC) :");
  348.     scanf("%lf",&hour);
  349.  
  350.  /*  hour = hour + 5;    make up difference for UTC time 
  351.     if ( hour > 24 )  
  352.      hour = hour - 24;     can't have 25 hours corection for over 24 
  353.  */
  354.    StartTime += hour/24;
  355.  
  356.     printf("Duration (Days) :");
  357.     scanf("%lf",&duration);
  358.     EndTime = StartTime + duration;
  359.  
  360.     printf("Time Step (Minutes) :");
  361.     scanf("%lf",&StepTime);
  362.     StepTime /= MinutesPerDay;
  363. }
  364.  
  365. PrintMode(OutFile,Phase)
  366. FILE *OutFile;
  367. {
  368.     int CurMode;
  369.  
  370.     for (CurMode = 0; CurMode < NumModes; CurMode++)
  371.     if ((Phase >= Modes[CurMode].MinPhase
  372.         && Phase <= Modes[CurMode].MaxPhase)
  373.           || ((Modes[CurMode].MinPhase > Modes[CurMode].MaxPhase)
  374.           && (Phase >= Modes[CurMode].MinPhase
  375.             || Phase <= Modes[CurMode].MaxPhase)))
  376.         {
  377.         fprintf(OutFile,"    %s",Modes[CurMode].ModeStr);
  378.         break;
  379.         }
  380. }
  381.  
  382.  
  383. main()
  384. {
  385.     double ReferenceOrbit;    /* Floating point orbit # at epoch */
  386.     double CurrentTime;     /* In Days               */
  387.     double CurrentOrbit;
  388.     double AverageMotion,CurrentMotion;
  389.     double MeanAnomaly,TrueAnomaly;
  390.     double SemiMajorAxis;
  391.     double Radius;        /* From geocenter           */
  392.     double RAANPrecession,PerigeePrecession;
  393.     double SSPLat,SSPLong;
  394.     long   OrbitNum,PrevOrbitNum;
  395.     int    Day,PrevDay;
  396.     double Azimuth,Elevation,Range;
  397.     double PrevRange,Velocity,Doppler;
  398.     int    Phase;
  399.     char   FileName[100];
  400.     FILE   *OutFile;
  401.     int    DidApogee;
  402.     double TmpTime,PrevTime;
  403.     int    PrevVisible;
  404.   
  405.     printf("\n");
  406.     printf("%s\n",Str1);
  407.     printf("%s\n",Str2);
  408.     printf("\n\n");
  409.  
  410.     GetSatelliteParams();
  411.     GetSiteParams();
  412.     GetSimulationParams();
  413.  
  414.     printf("Output file (Filename or * for CRT) :");
  415.     gets(FileName);    /* Skip previous RETURN */
  416.     gets(FileName);
  417.  
  418.  
  419.     if (strlen(FileName) > 0)
  420.     {
  421.     if ((OutFile = fopen(FileName,"w")) == 0)
  422.         {
  423.         printf("Can't write to %s\n",FileName);
  424.         exit(-1);
  425.         }
  426.     }
  427.       else OutFile = stdout;
  428.  
  429.     fprintf(OutFile,"%s Element Set %d\n",SatName,ElementSet);
  430.  
  431.     fprintf(OutFile,"%s\n",SiteName);
  432.  
  433.     fprintf(OutFile,"Doppler calculated for freq = %lf MHz\n",BeaconFreq);
  434.  
  435.     SemiMajorAxis = 331.25 * exp(2*log(MinutesPerDay/epochMeanMotion)/3);
  436.     GetPrecession(SemiMajorAxis,Eccentricity,Inclination,&RAANPrecession,
  437.             &PerigeePrecession);
  438.  
  439.     ReferenceOrbit = EpochMeanAnomaly/PI2 + EpochOrbitNum;
  440.  
  441.     PrevDay = -10000; PrevOrbitNum = -10000;
  442.     PrevTime = StartTime-2*StepTime;
  443.  
  444.     BeaconFreq *= 1E6;        /* Convert to Hz */
  445.  
  446.     DidApogee = 0;
  447.  
  448.     /* Start the loop one step early, to init OldRange */
  449.     for (CurrentTime = StartTime-StepTime; CurrentTime <= EndTime;
  450.         CurrentTime += StepTime)
  451.     {
  452.  
  453.     /* Fix these later to correct for drag    */
  454.     AverageMotion = epochMeanMotion;
  455.     CurrentMotion = epochMeanMotion;
  456.     SemiMajorAxis = 331.25 * exp(2*log(MinutesPerDay/CurrentMotion)/3);
  457.  
  458.     CurrentOrbit = ReferenceOrbit +
  459.             (CurrentTime-EpochDay)*AverageMotion;
  460.     OrbitNum = CurrentOrbit;
  461.  
  462.     MeanAnomaly = (CurrentOrbit-OrbitNum)*PI2;
  463.  
  464.     TmpTime = CurrentTime;
  465.     if (MeanAnomaly < PI)
  466.         DidApogee = 0;
  467.     if (PrintApogee && !DidApogee && MeanAnomaly > PI)
  468.         {            /* Calculate Apogee */
  469.         TmpTime -= StepTime;   /* So we pick up later where we left off */
  470.         MeanAnomaly = PI;
  471.         CurrentTime=EpochDay+(OrbitNum-ReferenceOrbit+0.5)/AverageMotion;
  472.         }
  473.  
  474.     TrueAnomaly = Kepler(MeanAnomaly,Eccentricity);
  475.     Radius = SemiMajorAxis*(1-SQR(Eccentricity))
  476.             / (1+Eccentricity*cos(TrueAnomaly));
  477.  
  478.     GetSubSatPoint(EpochDay,EpochRAAN,EpochArgPerigee,
  479.        Inclination,Eccentricity,RAANPrecession,PerigeePrecession,
  480.        CurrentTime,TrueAnomaly,&SSPLat,&SSPLong);
  481.     GetBearings(SiteLat,SiteLong,SiteAltitude,SSPLat,SSPLong,Radius,
  482.         &Azimuth,&Elevation,&Range);
  483.  
  484.     Velocity = (Range-PrevRange)/((CurrentTime-PrevTime)*SecondsPerDay);
  485.     PrevRange = Range;
  486.  
  487.     if (Elevation >= SiteMinElev && CurrentTime >= StartTime)
  488.         {
  489.         Day = CurrentTime + HalfSecond;
  490.         if (((double) Day) > CurrentTime+HalfSecond)
  491.         Day -= 1;    /* Correct for truncation of negative values */
  492.  
  493.         if (OrbitNum == PrevOrbitNum && Day == PrevDay && !PrevVisible)
  494.         fprintf(OutFile,"\n");    /* Dipped out of sight; print blank */
  495.  
  496.         if (OrbitNum != PrevOrbitNum || Day != PrevDay)
  497.         {            /* Print Header */
  498.                 fprintf(OutFile,"\n");
  499.                 PrintDate(OutFile,Day);
  500.                    fprintf(OutFile,"  ----Orbit # %ld-----\n",OrbitNum);
  501.               fprintf(OutFile," U.T.C.   Az   El  Doppler  Range");
  502.         fprintf(OutFile,"  Height  Lat  Long  Phase(%3.0lf)\n",
  503.                 MaxPhase);
  504.         }
  505.         PrevOrbitNum = OrbitNum; PrevDay = Day;
  506.             PrintTime(OutFile,CurrentTime + HalfSecond);
  507.  
  508.         Doppler = -BeaconFreq*Velocity/C;
  509.  
  510.         fprintf(OutFile,"  %3.0lf  %3.0lf  %6.0lf  %6.0lf",
  511.         Azimuth*DegreesPerRadian,
  512.         Elevation*DegreesPerRadian,Doppler,Range);
  513.  
  514.         fprintf(OutFile,"  %6.0lf  %3.0lf  %4.0lf",
  515.         Radius - EarthRadius,SSPLat*DegreesPerRadian,
  516.         SSPLong*DegreesPerRadian);
  517.  
  518.         Phase = (MeanAnomaly/PI2*MaxPhase + perigeePhase);
  519.         while (Phase < 0)
  520.         Phase += MaxPhase;
  521.         while (Phase >= MaxPhase)
  522.         Phase -= MaxPhase;
  523.  
  524.         fprintf(OutFile,"  %4d", Phase);
  525.         PrintMode(OutFile,Phase);
  526.         if (PrintApogee && (MeanAnomaly == PI))
  527.         fprintf(OutFile,"    Apogee");
  528.         fprintf(OutFile,"\n");
  529.         PrevVisible = 1;
  530.         }
  531.      else
  532.         PrevVisible = 0;    
  533.     if (PrintApogee && (MeanAnomaly == PI))
  534.         DidApogee = 1;
  535.     PrevTime = CurrentTime;
  536.     CurrentTime = TmpTime;
  537.     }
  538.     fclose(OutFile);
  539. }
  540.    END; 
  541.  
  542.  
  543.